# Execute R code to generate Figure 4 and Appendix Table Q1.
# Baron et. al (2024), Forthcoming, Political Behavior
# Couples Therapy for a Divided America: 
# Assessing the Effects of Reciprocal Group Reflection # on Partisan Polarization

# Note: Must run after executing the following STATA do files.
# baseline_clean.do, midline_clean.do, endline_clean.do, BArct_dataprep.do
# BArct_datapre.do should generate the cleaned final data file required for this analysis:
# Full_BArct_Data_panel.dta


rm(list=ls()) 
library(foreign)
library(ggplot2)
library(readstata13)
library(estimatr)
library(stargazer)
library(lubridate)
library(readr)
library(stringr)
library(lfe)
library(tidyverse)
library(ggridges)
library(broom)
library(ggsci)
library(forcats)
library(sjPlot)
library(ggeffects)
library(sjmisc)
library(jtools)
library(lm.beta)
library(ggpubr)



if(Sys.info()['user']=='dannychoi'){pathData<-"/Users/dannychoi/Dropbox/Projects/Better Angels RCT/10_Analysis/03_Code/Replication"}
setwd(pathData)
getwd()



### 1. Loading data sets ####

fulldat <- read.dta13(file="Full_BArct_Data_panel.dta")


### 2.  Attrition Analysis (IPW) ####

table(fulldat$attrit_mid) # midline attrition
table(fulldat$attrit_end) # endline attrition

table(fulldat$attrit_mid, fulldat$attrit_end)


#########################################
## Appendix Table Q.1: Attrition Table ##
#########################################

fulldat$attrit_midline <- NA
fulldat$attrit_midline <- ifelse(fulldat$part_base==1 & fulldat$part_mid==0, 1, fulldat$attrit_midline)
fulldat$attrit_midline <- ifelse(fulldat$part_base==1 & fulldat$part_mid==1, 0, fulldat$attrit_midline)
table(fulldat$attrit_midline)


fulldat$attrit_endline <- NA
fulldat$attrit_endline <- ifelse(fulldat$part_base==1 & fulldat$part_end==0, 1, fulldat$attrit_endline)
fulldat$attrit_endline <- ifelse(fulldat$part_base==1 & fulldat$part_end==1, 0, fulldat$attrit_endline)
table(fulldat$attrit_endline)

subdat2<- fulldat %>%
  filter(wave==1) %>%
  dplyr::select(treatment, wave, block, attrit_midline, attrit_endline, uid,
                b_affpol_noneg_comf, b_iat, b_therm_diff_std, 
                b_outparty_against_std, b_trust_diff_std, 
                b_outparty_friends_std, b_outparty_neigh_std, b_outparty_marry_std, 
                b_m1a, b_m2a, b_m2b, b_m3, b_m4, b_m5, b_m6, b_m7, b_ideoextreme, 
                b_idx_pol_engage, b_idx_outparty_contact, b_idx_empathy)

subdat3<- fulldat %>%
  filter(wave==1) %>%
  dplyr::select(treatment, wave, block, attrit_midline, attrit_endline, uid,
                b_affpol_noneg_comf, b_iat, b_therm_diff_std, 
                b_outparty_against_std, b_trust_diff_std, 
                b_outparty_friends_std, b_outparty_neigh_std, b_outparty_marry_std, 
                b_m1a, b_m2a, b_m2b, b_m3, b_m4, b_m5, b_m6, b_m7, b_ideoextreme, 
                b_idx_pol_engage, b_idx_outparty_contact, b_idx_empathy)


midline_int <- felm(attrit_midline ~ 
                      treatment*b_affpol_noneg_comf +
                      treatment*b_iat +
                      treatment*b_therm_diff_std + 
                      treatment*b_outparty_against_std + 
                      treatment*b_trust_diff_std + 
                      treatment*b_outparty_friends_std + 
                      treatment*b_outparty_neigh_std + 
                      treatment*b_outparty_marry_std + 
                      treatment*b_m1a + 
                      treatment*b_m2a + 
                      treatment*b_m2b + 
                      treatment*b_m3 + 
                      treatment*b_m4 + 
                      treatment*b_m5 + 
                      treatment*b_m6 + 
                      treatment*b_m7 + 
                      treatment*b_ideoextreme + 
                      treatment*b_idx_pol_engage + 
                      treatment*b_idx_outparty_contact + 
                      treatment*b_idx_empathy | factor(block) | 0 | 0, data=subset(fulldat, wave==1))
summary(midline_int)

endline_int <- felm(attrit_endline ~ 
                      treatment*b_affpol_noneg_comf +
                      treatment*b_iat +
                      treatment*b_therm_diff_std + 
                      treatment*b_outparty_against_std + 
                      treatment*b_trust_diff_std + 
                      treatment*b_outparty_friends_std + 
                      treatment*b_outparty_neigh_std + 
                      treatment*b_outparty_marry_std + 
                      treatment*b_m1a + 
                      treatment*b_m2a + 
                      treatment*b_m2b + 
                      treatment*b_m3 + 
                      treatment*b_m4 + 
                      treatment*b_m5 + 
                      treatment*b_m6 + 
                      treatment*b_m7 + 
                      treatment*b_ideoextreme + 
                      treatment*b_idx_pol_engage + 
                      treatment*b_idx_outparty_contact + 
                      treatment*b_idx_empathy | factor(block) | 0 | 0, data=subset(fulldat, wave==1))
summary(endline_int)


stargazer(midline_int, endline_int,
          type="text", single.row=T, keep.stat=c("adj.rsq", "n", "f"),
          covariate.labels = c("Treatment", 
                               "Explicit Affective Polarization", 
                               "Implicit Affective Polarization",
                               "Partisan Feeling Thermometer Gap",
                               "Negative Partisanship", 
                               "Partisan Trust Gap",
                               "Comfort w/ Out Party Friends",
                               "Comfort w/ Out Party Neighbors",
                               "Comfort w/ Out Party Marriage",
                               "Outgroup Stereotyping",
                               "Ideological Constraint",
                               "Ideological Extremity",
                               "Perceived Mass Ideo. Polarization",
                               "Outgroup Empathy",
                               "Outgroup Humanization",
                               "Moralization of Politics",
                               "Ingroup Identity Salience",
                               "Ideological Self-Placement",
                               "Political Engagement", 
                               "Out Party Countact",
                               "Dispositional Empathy",
                               "Treatment x Explicit Affective Polarization", 
                               "Treatment x Implicit Affective Polarization",
                               "Treatment x Partisan Feeling Thermometer Gap",
                               "Treatment x Negative Partisanship", 
                               "Treatment x Partisan Trust Gap",
                               "Treatment x Comfort w/ Out Party Friends",
                               "Treatment x Comfort w/ Out Party Neighbors",
                               "Treatment x Comfort w/ Out Party Marriage",
                               "Treatment x Outgroup Stereotyping",
                               "Treatment x Ideological Constraint",
                               "Treatment x Ideological Extremity",
                               "Treatment x Perceived Mass Ideo. Polarization",
                               "Treatment x Outgroup Empathy",
                               "Treatment x Outgroup Humanization",
                               "Treatment x Moralization of Politics",
                               "Treatment x Ingroup Identity Salience",
                               "Treatment x Ideological Self-Placement",
                               "Treatment x Political Engagement", 
                               "Treatment x Out Party Countact",
                               "Treatment x Dispositional Empathy"), 
          out="attrition_interactions.tex")









# Logistic regression to estimate probabilities of being observed


subdat2$observed <- ifelse(subdat2$attrit_midline==1, 0, 1)

subdat2 <- subdat2 %>%
  filter(!is.na(b_affpol_noneg_comf)) %>%
  filter(!is.na(b_iat)) %>%
  filter(!is.na(b_outparty_against_std)) %>%
  filter(!is.na(b_outparty_neigh_std)) %>%
  filter(!is.na(treatment))

subdat2$pobs2 <- glm(observed ~ treatment*b_affpol_noneg_comf + treatment*b_iat + treatment*b_outparty_against_std + treatment*b_outparty_neigh_std, data=subdat2, family=binomial(link="logit"))$fitted

hist(subdat2$pobs2)

# Endline

subdat3$observed <- ifelse(subdat3$attrit_endline==1, 0, 1)

subdat3 <- subdat3 %>%
  filter(!is.na(b_affpol_noneg_comf)) %>%
  filter(!is.na(b_iat)) %>%
  filter(!is.na(b_therm_diff_std)) %>%
  filter(!is.na(b_m1a)) %>%
  filter(!is.na(b_m2a)) %>%
  filter(!is.na(b_m4)) %>%
  filter(!is.na(b_idx_empathy)) %>%
  filter(!is.na(treatment))

subdat3$pobs3 <- glm(observed ~ treatment*b_affpol_noneg_comf + treatment*b_iat + treatment*b_therm_diff_std +  
                       treatment*b_m1a + treatment*b_m2a +treatment*b_m4 + treatment*b_idx_empathy, data=subdat3, family=binomial(link="logit"))$fitted

hist(subdat3$pobs3)

# Generating weights
subdat2$wt <- 1/subdat2$pobs2
subdat3$wt <- 1/subdat3$pobs3

# Creating data set for midline IPW analysis

mid_wts <- subdat2 %>%
  dplyr::select(uid, pobs2, wt, observed)

end_wts <- subdat3 %>%
  dplyr::select(uid, pobs3, wt, observed)

mid_ipw <- inner_join(fulldat, mid_wts, by="uid")
mid_ipw <- mid_ipw %>%
  filter(wave==2)

end_ipw <- inner_join(fulldat, end_wts, by="uid")
end_ipw <- end_ipw %>%
  filter(wave==3)


# Restrict analysis to observed subjects
midaffpol <- lm(affpol_noneg_comf ~ treatment + factor(block), weights=wt, data=subset(mid_ipw, observed==1))
endaffpol <- lm(affpol_noneg_comf ~ treatment + factor(block), weights=wt, data=subset(end_ipw, observed==1))
midiat    <- lm(iat    ~ treatment + factor(block), weights=wt, data=subset(mid_ipw, observed==1))
endiat    <- lm(iat    ~ treatment + factor(block), weights=wt, data=subset(end_ipw, observed==1))
middonate <- lm(donate ~ treatment + factor(block), weights=wt, data=subset(mid_ipw, observed==1))
enddonate <- lm(donate ~ treatment + factor(block), weights=wt, data=subset(end_ipw, observed==1))

midaffpol  <- tidy(midaffpol)
endaffpol  <- tidy(endaffpol)
midiat     <- tidy(midiat   )
endiat     <- tidy(endiat   )
middonate  <- tidy(middonate)
enddonate  <- tidy(enddonate)

beta <-         midaffpol$estimate[2]
beta <- c(beta, endaffpol$estimate[2])
beta <- c(beta, midiat$estimate[2])
beta <- c(beta, endiat$estimate[2])
beta <- c(beta, middonate$estimate[2])
beta <- c(beta, enddonate$estimate[2])

se <-       midaffpol$std.error[2]
se <- c(se, endaffpol$std.error[2])
se <- c(se, midiat$std.error[2])
se <- c(se, endiat$std.error[2])
se <- c(se, middonate$std.error[2])
se <- c(se, enddonate$std.error[2])


fxplot <- cbind(beta, se)
fxplot <- as.data.frame(fxplot)

modelname <- c("Explicit", 
               "Explicit", 
               "Implicit", 
               "Implicit", 
               "Behavioral", 
               "Behavioral")

submodel <- c("Midline",
              "Endline",
              "Midline",
              "Endline",
              "Midline",
              "Endline")

fxplot <- cbind(fxplot, modelname, submodel)
fxplot$rb <- sprintf("%0.3f", round(fxplot$beta, digits = 4))


fxplot$modelname <- factor(modelname, 
                           level=c("Explicit", 
                                   "Implicit", 
                                   "Behavioral"))

# Separate plots for midline and endline
interval1 <- -qnorm((1-0.90)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)
interval3 <- -qnorm((1-0.99)/2)


fxplot_mid <- fxplot %>%
  filter(submodel=="Midline")

treatplot_mid <- ggplot(data=fxplot_mid, aes(x=fct_rev(modelname), y=beta, shape=factor(modelname))) +
  geom_pointrange(aes(y = beta, ymin = beta - se*interval2,
                      ymax = beta + se*interval2, linetype=factor(modelname)),
                  lwd = 1/3, position = position_dodge(width = 1/2), size=0.3) +
  xlab("") +
  ylab("Treatment effect") +
  theme_classic() +
  ylim(-2.0, 1.0) +
  geom_hline(yintercept=0, linetype="longdash", color = "black", linewidth=0.2) +
  theme(axis.text.x= element_text(size=10),
        text = element_text(family = "Times New Roman"),
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank(),
        axis.line=element_line(size=0.2),
        legend.position="bottom",
        plot.title = element_text(hjust = 0.5),
        legend.text=element_text(size=12)) +
  coord_flip() +
  scale_shape_manual(name=NULL, values=c(15, 17, 16)) +
  scale_linetype_manual(name=NULL,
                        values = c("solid", "longdash", "dashed")) +
  labs(title = "Midline") + 
  guides(shape = guide_legend(override.aes = list(linetype = "blank", size=0.5))) 
print(treatplot_mid)


# Separate plots for midline and endline
fxplot_end <- fxplot %>%
  filter(submodel=="Endline")

treatplot_end <- ggplot(data=fxplot_end, aes(x=fct_rev(modelname), y=beta, shape=factor(modelname))) +
  geom_pointrange(aes(y = beta, ymin = beta - se*interval2,
                      ymax = beta + se*interval2, linetype=factor(modelname)),
                  lwd = 1/3, position = position_dodge(width = 1/2), size=0.3) +
  xlab("") +
  ylab("Treatment effect") +
  theme_classic() +
  ylim(-2.0, 1.0) +
  geom_hline(yintercept=0, linetype="longdash", color = "black", linewidth=0.2) +
  theme(axis.text.x= element_text(size=10),
        text = element_text(family = "Times New Roman"),
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank(),
        axis.line=element_line(size=0.2),
        legend.position="bottom",
        plot.title = element_text(hjust = 0.5),
        legend.text=element_text(size=12)) +
  coord_flip() +
  scale_shape_manual(name=NULL, values=c(15, 17, 16)) +
  scale_linetype_manual(name=NULL,
                        values = c("solid", "longdash", "dashed")) + 
  labs(title = "Endline") + 
  guides(shape = guide_legend(override.aes = list(linetype = "blank", size=0.5))) 
print(treatplot_end)

# Combine for final figure 
library(patchwork)
combined <- treatplot_mid + treatplot_end & theme(legend.position = "bottom")
combined <- combined + plot_layout(guides = "collect")
combined
ggsave('figures/attrit_IPW_combined.png', combined, width=10, height=5)


